Breakdown of the Stokes-Einstein Relation in Supercooled Water 
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ABSTRACT 

Supercooled water exhibits a breakdown of the Stokes-Einstein relation between 
the diffusion constant D and the alpha relaxation time r a . For water simulated 
with the TIP5P and ST2 potentials, we find that the temperature of the decou- 
pling of diffusion and alpha relaxation correlates with the temperature of the 
maximum in specific heat that corresponds to crossing the Widom line T\y(P). 
Specifically, we find that our results for Dr a /T collapse onto a single master 
curve if temperature is replaced by T — T\y{P), where T\y{P) is the temperature 
where the constant-pressure specific heat achieves a maximum. Also, we find 
agreement between our ST2 simulations and experimental values of Dr a /T. We 
further find that the size of the mobile molecule clusters (dynamical hetero- 
geneities) increases sharply near T\y(P). Moreover, our calculations of mobile 
particle cluster size (n(t*)) w for different pressures, where t* is the time for which 
the mobile particle cluster size is largest, also collapse onto a single master curve 
if T is replaced by T — T W (P). The crossover to a more locally structured low den- 
sity liquid (LDL) environment as T — ► Tw(P) appears to be well correlated with 
both the breakdown of the Stokes-Einstein relation and the growth of dynamic 
heterogeneities. 

A 17th century study of the density maximum at 4°C [l| demonstrates the long history 
of water science . Since that time, dozens of additional anomalies of water have 
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been discovered jfj, including the sharp increase upon cooling of both the constant-pressure 
specific heat Cp and the isothermal compressibility Kp. These anomalies of water become 
more pronounced as water is supercooled. To explain these properties, a liquid-liquid (LL) 
critical point has been proposed [7]. Emanating from this critical point there must be loci of 
extrema of thermodynamic response functions such as Cp and Kp- These loci must coincide 
as the critical point is approached, since response functions are proportional to powers of 
the correlation length, and the locus of the correlation length maxima asymptotically close 
to the critical point defines the Widom line. 

In the supercooled region of the pressure-temperature phase diagram, the dynamic prop- 



erties of water show dramatic changes [8|, . 
is the Stokes-Einstein (SE) relation 



One basic relation among dynamic properties 
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where D is the diffusion constant, T is the temperature, kp is the Boltzmann constant, rj 
is the viscosity and a is the effective hydrodynamic radius of a molecule. This expression 
provides a relation between mass and momentum transport of a spherical object in a viscous 
medium. The SE relation describes nearly all fluids at T > 1.2 — 1.6 T s , where T g is the glass 
transition temperature, and since the hydrodynamic radius a is roughly constant, Drj/T is 



approximately independent of T 
Drj/T is no longer a constant 
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however, in most liquids, for T < 1.2 — 1.6 T, 
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water, the breakdown of the SE relation occurs about 1.8T g , in the same temperature region 



22| . For the case of 
3 temper, 



in which many of the unusual thermodynamic features of water occur [5|, [8|, |23j. 

Our aim is to evaluate to what degree the SE breakdown can be correlated with the 
presence of thermodynamic anomalies and the onset of spatially hetero gene ous dynamics, 



and how these features relate to the location of the Widom line |16l . [24j, |25|, |26|, |27| . From 
prior studies of water, we can already form an expectation for the correlation between the 
SE breakdown and the Widom line by combining three elements: (i) the Widom line is 
approximately known from the extrapolated power-law divergence of Kp 28J; (ii) the locus 
of points Tp)(P) where D extrapolates to zero is also known, and nearly coincides with 
Tw(P) at low pressures (see Fig.l of Ref. [29(); (iii_) the SE relation has been found to 
fail in liquids generally at the temperature Tp>(P) 30j. Combining these three results, one 



might not be surprised if the breakdown of the SE relation should occur near to the Widom 
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line for P < Pc, and it should continue to follow Tb(-P) for P > Pc- We will see that 
our results are consistent with this expectation, but reveal some unexpected insights. To 
address these issues, we perform molecular dynamics (MP) simulations of N = 512 waterlike 
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32|. which exhibits a LL critical point 



33[| . We carry out simulations in 



molecules interacting via the TIP5P potential 
at approximately T c « 217 K and P c « 340 MPa |32. 
the NPT ensemble for three different pressures P = MPa, 100 MPa, and 200 MPa, and 
for temperatures T ranging from 320 K down to 230 K for P = and 100 MPa, and down 
to 220 K for 200 MPa. We also anal yze MD simulations of N = 1728 waterlike molecules 
interacting via the ST2 potential 34|, [35] , which displays a LL critical point at Tc ~ 245 K 



and P = 180 MPa 



ensemble 



361 ] . The simulations for the ST2 model are performed in the NVT 



351 ] . To locate the Widom line, we first analyze the isobaric heat capacity Cp for 



the TIP5P model [Fig. [11(a)] . For the ST2 model, we use the results for the Widom line 



from Ref. 



36|. 



We next explore the possible relation between the Widom line and the breakdown of the 
SE relation(Eq. Cp). We first calculate the diffusion constant via its asymptotic relation to 
the mean-squared displacement, 
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lim 

t— >oo 



(Ir^-r^O)! 2 ) 
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where Tj(t) is the position of oxygen j at time t. It is difficult to accurately calculate the 
viscosity rj in simulations, so we instead calculate the alpha relaxation time r a (which is 
expected to have nearly the same T dependence as rf) as the time at which the coherent 
intermediate scattering function 

(p(q,t)p(-q,0)> 



F(q,t) 



<p(q,0)p(-q,0)>" 



(3) 



decays by a factor of e. Here p(q, t) = exp[—iq ■ TCj(t)] is the Fourier transform of the 
density at time t, Tj{t) is the position of oxygen j at time t, q is a wave vector and the 
brackets denote an average over all |q| = q and many equilibrium starting configurations. 
We calculate F(q, t) at the value of q of the first maximum in the static structure factor 
F(q,0). It is important that we use the coherent scattering function (as opposed to the 
incoherent, or self-scattering function), since we want to capture collective relaxation that 
contributes to 77. Hence the SE relation can be written as 

Dt, 



T 



constant. 



(4) 



We see that both r a and D (Fig. H(b)) display rapid changes at low T. 

Figures 12(a) and [21(b) show Dr a /T as a function of T. At high T, Dr a /T remains 
approximately constant. At low T, Dr a /T increases indicating a breakdown of the SE 
relation (jTJ), which occurs in the same T region where the thermodynamic anomalies develop, 
and near the Widom line 7V(P). To test if there is a correlation between the SE breakdown 
and Tw(P), we plot Dr a /T against T — Tw(P) [Figs. [2(c) and (d)] and find the unexpected 
result that all the curves for different pressures overlap within the limits of our accuracy for 
both TIP5P and ST2 potentials. Hence Dr a is a function only of T — 2V(P), from which it 
follows that the locus of the temperature of the breakdown of the SE relation is correlated 
with T W (P). 

To better quantify the temperature where the SE relation breaks down, we use the fact 
that when the SE relation fails it can be replaced by a "fractional" SE relation D(r a /T)^ = 



const 
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351 ] . By plotting parametrically log(D) as a function of log(r Q /T), one can identify 
the crossover temperature T X (P) between the two regions by the intersection of the high 
T SE behavior and the low T fractional SE behavior [Fig. [3]. We confirm that the same 
collapse of Dr a /T can be found by replacing T with T — T x (P), demonstrating (Fig. [3]) that 
the loci of SE breakdown defined by T X (P) and Ty/(P) track each other for P < Pq. There 
is also some difference between the ST2 and TIP5P models in the relative location of the 
breakdown of the SE relation and T W (P), evidenced by the fact that the magnitude of the 
SE breakdown is different at T = 7V(P). As a result, the data for the two models will not 
collapse when plotted together. Moreover, we find that Dr]/T from the experimental data 
plotted against T — 7V(P) [Fig. [2(d)], coincides with the ST2 results, sug gest ing that the 
data collapse should exist for real water when T is replaced by T — T W (P) 37|, [38]. 

The above analysis has primarily focused on the behavior for P < Pc, where there is 
a Widom line. For P > Pq, water behaves similarly to simple glass forming liquids, and 
we expect the breakdown of the SE relation for water to be similar to other simple liquids. 
Specifically, the SE relation should break down at T « 1.2 — 1.6T S 13]. This breakdown 
approximately coincides with the temperature where the mode coupling description of the 



dynamics fails 
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30J. For the SPC/E model of water, the mode coupling temperature 



(Tmct(P)) locus has been evaluated, and the slope in the P-T plane of this locus changes 
from negative to positive for P > Pc [39] (the slope is positive for simple liquids). Corre- 
spondingly, there is also a decoupling of D and r a near Tmct{P) [39]. We find the same 
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behavior for the ST2 model. 

Since we find a correlation between Ty/{P) and the breakdown of the SE relation, the 
hypothesized connection between the SE breakdown and the onset of dynamical hetero- 
geneities (DH) suggests a connection between Tw(P) and the onset of DH. To investigate 
the behavior of the dynamic heterogeneities, we study the clusters formed by the 7% most 
mobile molecules [4^], defined as molecules with the largest displacements during a certain 
interval of time of length t. The clusters of the most mobile molecules are defined as follows. 
If in a pair of the most mobile molecules determined in the interval [to ,to + t], two oxygens 
at time to are separated by less than the distance corresponding to position of the first 
minimum in the pair correlation function (0.315 nm in TIP5P and 0.350 nm in ST2), this 
pair belongs to the same cluster. 

The weight averaged mean cluster size 

/ / ^ ( n2 (t)) 

measures the cluster size to which a randomly chosen molecule belongs, where (n(t)) is the 
number averaged mean cluster size. We show (n(t)) w in Fig. H]for TIP5P as a function of 
observation time interval t for different T at P = MPa [Fig. H](a)]. The behavior at higher 
P is qualitatively the same [Fig. IU(b)]. At low T, (n(t)) w has a maximum at the time t* 
associated with the breaking of the cage formed by the neighboring molecules (see 41[ and 
the references therein). Both the magnitude and the time scale t* of the peak grow as T 
decreases. At high T, this peak merges and becomes indistinguishable from a second peak 
with fixed characteristic time ~ 0.5 ps. By evaluating the vibrational density of states, we 
associate this feature with a low frequency vibrational motion of the system, probably the 



0-0-0 bending mode 



42|. 



To probe the temperature dependence of {n(t)) w , we plot the peak value {n(t*)) w in 
Fig. He) as a function of T for P = MPa, 100 MPa, and 200 MPa for the TIP5P model. 
At high T, (n(t*)) w is nearly constant, since at high T, clusters result simply from random 
motion of the molecules. Upon cooling near the Widom line, {n(t*)) w increases sharply. 
When (n(t*)) w is plotted as a function of T — 2V(P) [see Fig. Hfd)], we find that (similarly 
to the behavior of Dr a /T) the three curves for P = MPa, P = 100 MPa, and P = 200 MPa 
overlap, and it is apparent that the pronounced increase in (n(t*)) w occurs for T « T W (P). 

To further test that the breakdown of the SE relation in water is associated with the 
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onset of large DH near the Widom line, we show (n(t)) w for different T along an isochore 
of density p = 0.83 g/cm 3 for the ST2 model of water in Fig. 0(a). We show an isochore 
because only isochoric data are available from Ref. 36|. As in the case of TIP5P, we find the 
emergence of a second time scale larger than 0.5 ps in (n(t)) w near the crossing of the Widom 
line at this density. Similarly, {n(t*)) w increases sharply near the Widom line temperature 
[see Fig. Ep))]. Hence the sharp growth of DH and the appearance of a second time scale 
in (n(t)) w both occur near the Widom line. We also find that the magnitude of (n(t*)) w is 
larger for the ST2 model than for the TIP5P model at the Widom temperature, consistent 
with the above observation that the breakdown of the SE relation is more pronounced for 
the ST2 model than for the TIP5P model. 

Finally, we consider the correspondence between DH and static structural heterogeneity 
in the supercritical region; this region is characterized by large fluctuations spanning a wide 
range of structures, from HDL-like to LDL-like. To quantify these structural fluctuations, 
we calculate for the TIP5P model the temperature derivative of a local tetrahedral order 
parameter Q 43J. In Fig. Elfa), we show \(dQ/dT)p\ at different T for P = MPa, and 



100 MPa, and find maxima 



44) at 7V(P) |45|, |4fj|. The maxima in \{dQ/dT) P \ indicates 



that the change in local tetrahedrality is maximal at 2V(P), which should occur when 
the structural fluctuation of LDL-like and HDL-like components is largest. We see that 
the growth of the dynamic heterogeneity coincides with the maximum in fluctuation of 
the local environment. Also, since Q quantifies the orientational order, the fact that we 
find that \{dQ/dT)p\ has maximum at approximately the same temperature where Cp = 
T(dS/dT)p has a maximum, supports the idea that water anomalies are closely related to 
the orientational order present in water. 

Before concluding, we note that more than a dozen other phenomena have been correlated 
with Ty/(P). Some of these phenomena are from simulations, such as the crossover in the 
relaxation time of the fluctuations in orientational order parameter Q 47|] or the maximum 
in the temperature derivative of the number of hydrogen bonds per molecule 48|. Others 
are from only experiments, such as the sharp drop in the the temperature derivative of the 
zero-frequency structure factor observed by QENS or the appearance of a Boson peak, both 
observed by quasi-elastic neutron scattering [49|. Finally, some anomalies that correlate 
with the Widom line are found in both simulations and experiment, such as the dynamic 
"fragile-to-strong" crossover in the diffusion constant, or the sharp drop in temperature 
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derivative of of the mean squared displacement (seen in QENS, NMR, and simulations) 



44|, |45|, 15( 

In conclusion, we find that the breakdown of the Stokes-Einstein relation for P < Pq 
can be correlated with the Widom line. In particular, rescaling T by T — T\y(P) yields 
good data collapse of Dr a /T for different pressures. Rapid structural changes occur at T 
near the Widom line, where larger LDL "patches" emerge upon cooling. The size of the 
dynamic heterogeneities also increases sharply as the Widom line is crossed. The breakdown 
of the SE relation can be understood by the fact that diffusion at low T is dominated by 
regions of fastest moving molecules (DH) while the relaxation of the system as a whole is 
dominated by slowest moving molecules. Consistent with this, we found that the growth 
of mobile particle clusters occurs near the Widom line and also near the breakdown of the 

ius the SE breakdown in water is consistent with the LL-critical 
7J. Our results are also consistent with recent experimental 



SE relation for P < Pc- T 
point hypothesis 0, y, 4 
findings in confined water [8 
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FIG. 1: Temperature dependence of (a) Cp and (b) r a and D for P = MPa, 100 MPa and 
200 MPa. We note that while the temperature of the maximum decreases with pressure, the value 
of Cp(T) at the maximum also decreases. This could be an artifact of the TIP5P potential which 
is known not to correctly reproduce the radial distribution function at high pressures. 
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FIG. 2: (a) Dr a /T as a function of T for P = MPa, 100 MPa and 200 MPa for the TIP5P 
model. For all panels, Dr a /T is scaled by its high T value to facilitate comparison of the different 
systems, (b) Analog of Fig. [2ja) for the ST2 model, (c) Dr a /T as a function of T — T\y{P) for 
TIP5P. (d) Analog of Fig. [2] (c) for the ST2 model including the experimental curve of Dr//T for 
water. 
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FIG. 3: Locus in P-T plane of Cp ax and T x (P) for the ST2 model. Filled circle denotes the 
liquid-liquid critical point (Tc,Pc) in ST2 model. 
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FIG. 4: Mobile particle cluster size (n(t)) w (a) for different T at P = MPa and (b) for different 
P at T = 230 K. (c) (n(t*)) w for P = MPa, 100 MPa and 200 MPa. (d) (n(i*)) as a function of 
T — Tyy(P) for three different P. All these plots are for the TIP5P model. 
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FIG. 5: Dynamical heterogeneities in the ST2 model of water, (a) (n(t)) w as a function of t at 5 K 
intervals for p = 0.83 g/cm 3 from 265 K to 300 K. The bold line shows the T = 280 K isotherm 
where the constant volume specific heat Cy has a maximum, (b) (n(t*)) w as a function of T for 
p = 0.83 g/cm 3 . We indicate the temperature at which Cy has a maximum by a vertical arrow. 
The maximum of Cy occurs in the same region where other response functions, such as Cp, have 
maxima. 
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FIG. 6: \{dQ/dT) P \ as a function of T for P = MPa and 100 MPa and 200 MPa. We note 
that while the temperature of the maximum decreases with pressure, the value \{dQ/dT)p\ at 
the maximum for P = 200 MPa than the lower pressures. This could be an artifact of the 
TIP5P potential which is known not to correctly reproduce the radial distribution function at high 
pressures. 
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